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I .  INTRODUCTION 


In  free  flight  studies,  the  main  task  is  to  determine  from  obser¬ 
vations  of  a  missile's  motion  the  values  of  the  aerodynamic  coefficients 
appearing  in  the  differential  equations  describing  that  motion.  For 
normal  enclosed-range  firing  conditions  (symmetric  shell,  nearly  hori¬ 
zontal  flight,  small  yaw,  constant  or  only  slowly  varying  spin,  etc.), 
we  can  approximate  the  solution  to  the  differential  equations  quite 
adequately  by  convenient  closed-form  expressions.  These  expressions 
involve  certain  constants  directly  related  to  the  aerodynamic  coeffi¬ 
cients.  The  values  of  these  constants  are  determined  by  a  least  squares 
fit  of  the  observed  data  to  the  closed-form  expressions. 

For  the  past  two  decades,  the  general  technique  of  fitting  observed 
data  to  convenient  closed-form  expressions  has  been  the  heart  of  free- 
flight,  enclosed-range  data  reduction1*2  .  Highly  successful  results 
have  been  obtained  for  a  variety  of  missile  shapes  and  sizes.  Over  the 
Years,  the  technique  has  been  gradually  refined  and  extended  to  cover 
many  types  of  force  and  moment  nonlinearities.  Unfortunately,  such 
extension  often  requires  that  a  number  of  rounds  be  fired  at  different 
Mach  numbers,  yaw  levels,  and  so  on,  to  obtain  a  single  set  of  coeffi¬ 
cients.  The  process  can  be  costly  in  dollars  and  time.  Moreover,  the 
technique  can  be  stretched  only  so  far;  an  occasional  round  has  defied 
analysis  by  conventional  procedures. 

It  is  relatively  easy  in  these  troublesome  nonlinear  situations  to 
specify  -  by  experience  and  by  cunning  -  the  sort  of  nonlinear  terms 
that  must  appear  in  the  differential  equations  of  motion  in  order  to 
produce  the  observed  behavior.  It  is  much  more  difficult  to  find 
convenient  pseudo-solutions  that  will  (1)  represent  the  motion  adequately 
and  (2)  contain  constants  that  can  be  easily  related  to  the  coefficients 
of  the  differential  equations.  What  we  need  in  these  situations  is  a 
method  that  doesn't  require  any  knowledge  or  assumptions  on  our  part 
regarding  the  form  of  the  solution  to  the  differential  equation. 

The  problem  is  essentially  one  of  parameter  optimization.  We  are 

given 


a.  the  form  of  a  set  of  differential  equations  involving  unknown 
constant  parameters  (coefficients  and  initial  conditions); 

b.  a  set  of  first  estimates  for  the  parameters; 

c.  a  set  of  discrete  measurements  on  one  or  more  of  the  dependent 
variables ; 

d.  a  criterion  function  of  the  parameters  (sr.y,  the  sum  of  the 
squares  of  the  residuals  for  a  given  fit) . 


* References  are  listed  on  page  59. 
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The  problem  is  to  devise  a  routine  for  adjusting  the  parameters  auto- 
matically  so  as  to  minimize  the  criterion  function. 

The  problem  has  received  much  attention,  particularly  in  recent 
years  with  the  proliferation  of  high-speed  computers,  and  many  ingenious 
schemes  for  optimizing  the  parameters  have  been  proposed.  For  continuous 
rather  than  discrete  input  data,  Meissinger3 >**  approached  the  criterion 
function  minimum  by  a  path  of  approximately  steepest  descent,  using  an 
analog  computer.  On  the  other  hand,  the  techniques  devised  by  Goodman5"7 
and  by  the  team  of  Chapman  and  Kirk®  to  handle  discrete  data  are  better 
suited  to  the  digital  computer.  In  this  paper,  we  will  be  concerned 
primarily  with  the  Chapman-Kirk  technique.  The  Goodman  and  Chapman- 
Kirk  methods  are  quite  similar,  although  Chapman  and  Kirk  were  unaware 
of  Goodman's  work  when  they  presented  their  own  results  at  the  AIM 
Seventh  Aerospace  Sciences  Meeting,  New  York,  1969.  From  the  stand¬ 
point  of  the  aerodynamicist.  Chapman  and  Kirk's  contribution  was  to 
apply  the  process  successfully  to  representative  aerodynamic  cases  and, 
perhaps  more  important,  to  present  their  results  at  the  meeting  and 
later  in  a  journal8  where  it  came  to  aerodynami cists '  attention.  (It 
is  unfortunate  but  often  true  that  when  a  pertinent  article  such  as 
Goodman's  appears  in  a  mathematical  journal,  the  aerodynamicist  either 
overlooks  it  or  fails  to  recognize  its  applicability  to  his  work.) 

Applications9'13  of  the  Chapman-Kirk  technique  are  growing  more 
and  more  sophisticated.  The  present  report  documents  one  such  appli¬ 
cation  carried  out  by  the  Armament  Department,  General  Electric  Company, 
Burlington,  Vermont,  for  the  U.  S.  Army  Aberdeen  Research  and  Development 
Center  (ARDC) ,  Aberdeen  Proving  Ground,  Maryland,  under  Government 
Contract  No.  DMD05-71-C-0265  during  the  period  4  February  to  20  April 
1971.  The  results  of  this  work  have  also  been  issued  as  a  General 
Electric  report1**. 

The  raw  data  for  this  study  consisted  of  time,  position  and  orien¬ 
tation  measurements  at  discrete  points  along  the  trajectory  for  twelve 
rounds  fired  in  the  Transonic  Free  Flight  Range15,  Ballistic  Research 
Laboratories  (BRLj .  The  twelve  rounds  consisted  of  four  spin-stabilized 
projectile  types  (see  Table  I  and  Figure  1).  Each  of  the  rounds  had 
been  previously  reduced11*16  by  the  usual  range  technique  (with  varying 
degrees  of  success)  and  each  was  hand-picked  for  the  present  assignment. 

Sonn  i  the  twelve  rounds  could  be  fitted  by  assuming  relatively 
simple  force  and  moment  systems;  these  rounds  were  chosen  to  enable  the 
Chapman-Kirk  technique  to  get  a  foot  in  the  door.  Other  rounds  of  the 
twelve  were  oddities  that  had  already  annoyed  and  frustrated  a  team  of 
data  analysts;  these  rounds  were  chosen  to  give  the  Chapman-Kirk  tech¬ 
nique  a  good  work-out.  The  primary  purpose  of  the  present  study  was 
to  see  how  well  the  Chapman-Kirk  technique  could  determine  the  values 
of  the  nonlinear  coefficients  present  in  the  force  and  moment  expres¬ 
sions. 
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II.  THE  CHAPMAN-KIRK  TECHNIQUE 


Suppose  for  the  moment  that  we  have  a  system  defined  by  elementary 
equations*.  Let 

N1  =  the  number  of  measured  dependent  variables 

N23  =  the  number  of  unknown  constants  in  the  system 

N4  =  the  number  of  measurements  taken  on  each  of  the  N1  dependent 
variables 

where  N4  is  greater  than  N23.  (The  notation  here  is  not  entirely  capri¬ 
cious.  An  N2  and  N3,  lying  between  N1  and  N4,  will  be  introduced  later 
and  N23  will  be  the  sum  of  N2  and  N3.)  Assume  that  we  can  write  our 
system  of  equations  in  the  form 

Yi  “  U,  i  ^2*  *  *  ' 1  ^N2  3)  »  *  ”  1  ^ »  *  ..|N1  (1) 

where  y^  is  the  i-th  dependent  variable  (i  =  1,  2,  .  .  .  Nl) 
f^  is  a  single-valued  elementary  function 
x  is  The  independent  variable 

a^  is  the  n-th  unknown  constant  parameter  (n  =  1,  2,  .  .  .,  N23) 

We  are  given  a  set  of  measurements,  which  we  represent  by  the  Nl  x  NJ 
matrix  D  =  (d^m) ,  and  a  vector  X  =  (xm)  of  the  corresponding  N4  values 

of  the  independent  variable.  That  is, 

d.  =  the  measured  value  of  y.  at  x 
1m  '  1  m 

i  =  1,  2,  ....  Nl 

m  =  1 ,  2 ,  .  .  . ,  N4 

For  any  parameter  vector  t  =  (a^) ,  we  can  obtain  from  (1)  the  correspond¬ 
ing  solution  y.  (x  )  at  each  point  x  .  The  problem  is  to  determine  the 

value  of  A  that  minimizes  e,  the  sum  of  the  squares  of  the  residuals: 


*By  an  elementary  equation,  we  mean  an  equation  involving  only  elementary 
functions  and  a  finite  nunber  of  arithmetical  operations .  Specifically , 
we  are  excluding  differential  equations . 
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e 


(2) 


N4  N1 

■  E  E  [di-  •  (x-)]2 

m=l  i-1 


(For  convenience,  we  omit  weighting  factors  in  (2)  and  the  succeeding 
discussion;  such  factors  might  be  needed,  for  example,  to  insure  that 
the  terms  in  (2)  are  dimensionally  equal.)  Now  e  will  be  at  a  minimum 
only  when  its  partial  derivatives  with  respect  to  each  of  the  para¬ 
meters  is  zero.  We  introduce  a  convenient  notation  for  these  partial 
derivatives : 


3  yA 

p.  =  — - 
in  9  a 

n 

i  -  1,  2 . N1 

n  =  1,  2,  .  .  .,  N23 


(3) 


Because  they  reflect  the  influence  of  each  parameter  change,  the  Pip 

are  sometimes  called  "influence"  (or  "sensitivity")  coefficients.  For 
a  minimum  e,  we  must  have 


9  e 

rq  ' 


N4  N1 

2  E  E  Idi.-yi  (X"H  Pin-° 

m=  1  i=  1 


(4) 


Equations  (4)  constitute  a  set  of  N23  equations  in  the  N23  unknown  param¬ 
eters;  the  values  of  the  parameters  satisfying  (4)  are  the  desired 
optimum  values. 

If  the  functions  f^  are  linear  in  the  parameters  a^: 


f. 

l 


a  (x)  +  a  <J>.  (x)  +  .  . 

l  xl  2  12 


*  ^23  *iN23  (x) 


then 


P. 

in 


rin 


(x) 


and  set  (4)  is  also  linear  in  the  parameters;  hence  (4)  is  easily 
solvable  (in  theory).  If  the  functions  f  are  nonlinear  in  the  para¬ 
meters.  owever,  then  solving  (4)  can  be  quite  difficult.  The  usual 
way  of  £  this  difficulty  is  to  approximate  the  variables  y^  by  their 
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linearly  truncated  Taylor  expansions  about  a  given  set  of  values  for 
the  parameters : 

N23 

yi *  h +  JZ  Kk  Aak 

k=l 


where  the  circumflex  (*)  denotes  evaluation  at  the  given  parameter 
values  and  where  Aa^  is  an  indicated  change  in  the  given  value  of  a^. 


If  (5)  is  substituted  in  (4) ,  and  if  the  in  (4)  are  replaced 

by  P.  ,  we  have 
1  in’ 

N4  N1 

m=l  i«l 

or 

N23 

7~.  «„k  4*k  '  »»■”■>• 2 . »23  <6> 

k=l 


N23 


k=l 


where 


nk 


P.  (x  )  P..  (x  ) 
in  m'  lk  m' 


(7) 


N4  N1 


-cr; 


m= 1  i»  1 


P. 

in 


(8) 


Equation  (6)  represents  a  set  of  N23  linear  equations  in  the  ^23  incre¬ 
ments  Aan>  Hence  we  can  solve  (6)  for  the  increments: 

N23 

4an  *  YZ  V  "  ■  *■  2 . N23  <9> 

k=l 


19 


where 


nk 


the  (n,  k)-th  element  of  the  inverse  of  matrix  (a^) 


cofactor  of  element  a  , 
_ _ _  nk 

determinant  of1  matrix  (o>nk) 


The  new  value  of  a  is  then  obtained  by  adding  Aa  to  the  old  value, 
n  n 


If  the  initial  estimates  of  the  parameters  are  "close  enough"  to 
the  "true"  va'ues,  (5)  will  be  an  adequate  approximation  and  the  new 
values  of  the  parameters  will  produce  a  smaller  criterion  function  c. 
The  process  can  then  be  repeated  as  many  times  as  necessary  until  some 
convergence  criterion  is  satisfied.  The  probable  error  of  the  fit  at 
the  end  of  any  iteration  is  given  by 


PE  =  0.6745 


r _ - _ r 

I  N1  x  N4  -  N23 


and  an  estimate  of  the  probable  error  in  a^  is  given  by 


(10) 


H  (an) 


(a  )**  •  pe 
iur 


(ID 


This,  in  brief,  is  the  well-known  process  of  "differential  corrections," 
as  applied  to  elementary  equations  that  are  nonlinear  in  the  parameters. 

For  a  system  of  ordinary  differential  equations,  the  situation  is 
naturally  more  complicated.  Assume  that  the  given  system  of  differential 
equations  has  been  reduced  to  a  system  of  N2  first  order,  possibly  non¬ 
linear  eouations*: 


d/j 

=  fj  (x’Vy2' 


’  yN2 ’  W 


yj  (xo>  =  *N3  ♦  j 


1,  2,  .  .  .,  N2 


*It  is  not  neoe88ary  when  carrying  out  the  Chapman-Kirk  technique  to 
reduce  the  given  system  to  first  order  equations;  this  was  done  here 
merely  to  aid  the  exposition. 
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where 


N2  =  the  total  number  of  dependent  variables  in  the  system 

N3  =  the  number  of  unknown  parameters  appearing  explicitly 
in  the  differential  equations 

and  where  N1  and  N4  are  defined  as  before.  Note  that  we  have  juggled 
the  notation  so  that  while  only  N3  parameters  appear  explicitly  in  (12), 
there  are  still  N23  (=  N2  +  N3)  parameters  to  be  determined: 

a^,  a  ,  .  .  a^j  (the  N3  explicit  parameters) 

^3+1’  ^3+2’  ‘  '  *»  ^23  N2  un*cnown  initial  conditions) 


where 


1  <  N1  <S  N2  <  N23  <  N4  (13) 

The  same  differential  correction  technique  that  was  applied  to  (1) 
can  be  applied  to  (12).  Equations  (6-11),  which  depend  only  on  the 
definition  (2)  of  the  criterion  function  e,  are  still  valid.  However, 
a  new  difficulty  arises  when  the  given  set  of  equations  are  ordinary 
differential  equations:  how  do  we  evaluate  the  dependent  variables 
and  their  partial  derivatives  appearing  in  (7)  and  (8)?  For  the  case 
of  non-differential  equations,  this  was  no  problem;  we  were  presumably 
given  explicit  expressions  for  each  of  the  variables  and  could  easily 
write  down  expressions  for  the  required  partials.  For  a  given  set  of 
differential  equations,  however,  we  will  not,  in  general,  know  the  form 
of  the  solution.  Values  of  the  dependent  variables  can  be  obtained  by 
some  numerical  integration  scheme,  but  part  of  the  problem  remains:  how 
do  we  obtain  the  required  values  of  the  partial  derivatives? 

Chapman  and  Kirk  tried  various  schemes  for  evaluating  these  partial 
derivatives  and  finally  settled  on  the  method*  of  "parametric  different¬ 
iation."  This  method  consists  of  formally  differentiating  the  given  set 
(12)  of  differential  equations  with  respect  to  each  of  the  N23  constants 
to  be  determined.  We  have 


* The  method  is  not  new.  It  was  used ,  for  example ,  by  the  previously 
cited  Meis3inger  and  Goodman  (and  apparently  by  Knalder^7 ,  with  whose 
work  we  are  unfamiliar  but  who  is  referenced  in  References  10  and  12). 
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or,  assuming  that  the  y.  are  continuous  in  x  and  an>  so  that  the  order 
of  differentiation  can  be  interchanged. 


d  P.  3  f. 

—liL  x  _ 1 

dx  3  a 

n 

j  =  1,  2,  .  .  .,  N2 

(14a) 

n  =  1 ,  2,  .  .  . ,  N23 

where 

Pjn  "  1 

(if  n  -  j  =  N3) 

=  0 

(otherwise) 

(14b) 

The  somewhat  strange-looking  initial  conditions  (14b)  merely  reflect  the 
fact  that  Pjn  is  initially  1  if  and  only  if  represents  the  initial 

value  of  yy  By  (12),  this  occurs  if  and  only  if  a^  ■  a^  +  j;  hence 

if  and  only  if  n  =  N3  +  j  . 

What  we  have  done  above  is  to  derive  an  auxiliary  set  of  N2  x  N3 
equations  (14)  whose  solutions  are  the  partials  P^n,  some  of  which 

(those  for  j  <  Nl)  are  needed  in  solving  (6).  For  a  given  set  of 
estimates  of  the  parameters  and  initial  conditions,  we  can  integrate  by 
some  numerical  scheme  both  the  original  set  of  equations  (12)  and  the 
auxiliary  set  (14).  (The  original  set  may  or  may  not  be  linear;  the 
auxiliary  set  will  always  be  linear.)  The  numerical  integration  yields 
the  values  of  the  dependent  variables  and  the  influence  coefficients 
required  to  solve  (6)  for  the  parameter  changes. 

Except  for  this  more  laborious  way  of  determining  y^  (xffl)  and 

Pin  (xm) ,  the  procedure  for  determining  the  unknown  constants  of  a 

set  of  differential  equations  is  the  same  as  for  non-differential 
equations. 

As  a  final  aside,  we  note  a  possible  future  improvement.  In  the 
N23-dimensional  parameter  space,  the  truncated  Taylor  series  technique 
proceeds  from  a  given  point  (whose  coordinates  are  the  given  estimates 
of  the  N23  constants)  in  the  direction  of  the  vector  AA  =  (Aa^) 

obtained  by  solving  (6).  The  method  of  steepest  descent,  on  the  other 
hand,  proceeds  from  the  given  point  in  the  direction  of  the  negative 
gradient  of  e.  Marquardt 1 8 » 1 9  points  out  that  these  two  directions 
are  nearly  perpendicular,  while  the  optimum  direction  lies  somewhere 
in  between.  To  proceed  in  approximately  the  optimum  direction,  he 
suggests  replacing  in  (6)  by 
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(15) 


'kn 


(1  +  X)  a^n  for  k  =  n 


a.  for  k  ¥■■  n 
kn 


where  X  is  a  fudge  factor  constant  whose  value  should  be  changed  from 
one  iteration  to  the  next  according  to  a  few  simple  rules  (the  rules 
are  listed  -  in  slightly  different  form  -  in  References  18  and  19). 
Using  Marquardt's  magic  X,  the  Chapman-Kirk  process  often  converges 
for  initial  guesses  far  outside  the  previous  region  of  convergence. 
The  Marquard  alg  rithm  was  pointed  out  to  us  by  Chapman  himself,  who 
has  used  it  (subsequent  to  the  work  reported  on  in  Reference  8)  with 
great  success  in  hitherto  intractable  cases.  The  algorithm  was  not 
used  in  the  investigation  covered  by  this  report  because  we  were  not 
aware  of  it  at  the  time. 


III.  INPUT  DATA,  COORDINATE  SYSTEMS  AND  YAW  VARIABLES 

A  missile  fired  in  the  BRL  Transonic  Range  is  observed  at  twenty- 
five  spark-photography  stations  distributed  along  a  680-foot  portion 
of  the  trajectory.  The  observed  data  at  each  station  consists  of 

a.  the  elapsed  time  t,  reckoned  from  the  instant  the  spark  at  the 
first  station  was  triggered  by  the  passing  missile.  The  time  error  in 
a  properly  functioning  timer  is  estimated  to  be  no  more  than  one  micro¬ 
second.  Only  about  two-thirds  of  the  stations  are  instrumented  at 
present  to  furnish  timing  data. 

b.  (x,y,z)R:  the  position  vector  of  the  missile's  CG  in  a  range 

(earth-fixed)  coordinate  system  XYZ.  The  X-axis  is  directed  down-range 
along  the  intersection  of  a  horizontal  plane  with  the  vertical  plane 
containing  the  gun;  the  Y-axis  lies  in  this  horizontal  plane,  directed 
to  the  left  of  an  observer  facing  downrange;  the  Z-axis  is  directed  up¬ 
ward.  The  error  in  any  position  measurement  should  be  no  greater  than 
0.003  meter. 


c.  (0,  C  .  C  )cn:  the  yaw  vector  in  a  fixed-plane  coordinate 

2  3 

system  XXX.  The  X  -axis  lies  along  the  missile's  longitudinal  axis, 
12  3  1 

directed  from  the  CG  (the  origin  of  the  fixed-plane  system)  to  the  nose 
the  X^-axis  is  constrained  to  lie  in  the  horizontal  plane,  directed  to 

the  left  of  an  observer  facing  in  the  direction  of  the  positive  X  -axis 

l 

the  X^-axis  is  directed  according  to  the  right-hand  rule  (that  is,  up¬ 
ward).  If  u^,u  ,u^  are  the  velocity  components  in  the  fixed-plane 
system,  then 
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u 

-  _1 
V 


(16) 


e 

2 


U 


i  » 
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where  V  is  the  magnitude  of  the  velocity  vector.  The 
(16)  appear  because  in  range  studies  the  yaw  angle  is 
the  velocity  vector  to  the  X  -axis*.  The  magnitude  6 

given  by 


minus  signs  in 
measured  from 
of  the  yaw  is 


-  sin  (Xj.  (17) 


where  cXj,  is  the  yaw  angle,  the  so-called  total  angle  of  attack.  The 

angular  measurements  that  yield  the  yaw  components  are  usually  accurate 
to  within  0.002  radian. 


Although  differential  equations  describing  the  yawing  motion  can  be 

written  in  terms  of  £  and  £  ,  we  found  it  more  convenient  to  work  with 

2  3 

the  related  Euler  angles  and  8.  These  angles  appear  in  the  trans¬ 
formation  matrix  that  converts  from  range  to  fixed-plane  coordinates . 

To  derive  this  matrix,  assume  the  existence  of  a  vector  3  extending 
from  the  range  system  origin  in  the  direction  of  the  positive  X^-axis 

(see  Figure  2a).  Rotate  the  range  system  XYZ  about  the  Z-axis  by  the 
angle  ^  so  that  X'  (the  rotated  X-axis)  coincides  with  the  projection 
of  the  vector  B  on  the  XY -plane.  The  magnitude  of  is  not  to  exceed 
180°;  if  this  requires  a  counterclockwise  rotation  about  the  Z-axis  (as 
in  Figure  2a),  \p  is  considered  positive  and  if  clockwise,  then  iji  is 
considered  negative.  Then  we  have 


*In  References  1  and.  2,  a  fixed-plane  system  X^X*X^  is  used ,  where  the 

X7 -  and  X'  -axes  have  opposite  directions  to  the  X  -  and  X  -axes,  respec- 
^  3  2  3 

tively.  In  this  X^X^X^  system,  the  yaw  angle  is  measured  from  the  X^- 

axis  to  the  velocity  vector.  If  u  ,u7  ,uy  are  the  velocity  components 

end  0,5'  are  the  yaw  components  in  the  X  X' X'  system,  then 

2  3  12  3 


u' 

JL  = 
V 


u 

.  JL 
V 


and  similarly ,  That  is,  the  yaw  components  in  the  two  fixed- 

plane  systems  are  identical. 
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cos  ip  sin  i|i  0\  /  x\ 

-  sin  ^  cos  i|/  Oil  y  I 

0  0  1  /  \  2  /  r  (18) 

Next,  rotate  the  intermediate  system  X'Y'Z'  about  the  Y' -axis  so  that  X" 
(the  rotated  X-axis)  coincides  with  vector  B.  Let  0  denote  the  angle 
from  V  to  8,  where  |e|  <  90°.  If  the  rotation  is  counterclockwise 
about  the  Y/-axis,  6  is  considered  positive;  if  the  rotation  is  clock¬ 
wise  (as  in  Figure  2a),  0  is  negative.  Then  we  have 

x"\  /  cos  0  0  -  sin  0W  x7  \ 

y"  =  I  o  i  o  y7 

Z"/  FP  y  sin  0  0  cos  0  /  \  z' /  R/  (19) 


This  final  system  X"Y"Z"  has  the  orientation  of  the  fixed-plane  system 
X^X  X^.  Thus  the  transformation  matrix  from  range  to  fixed-plane  co¬ 
ordinates  is  the  product  of  the  e -matrix  and  the  i^-matrix: 


COS  0  COS  1^ 

-  sin  [p 
sin  0  cos  tp 


cos  0  sin  tp 
cos  \p 
sin  0  sin 


(20) 


For  an  assumed  flat,  nonrotating  earth, 
has  the  form 


the  gravitational 


force  2 


5  =  (0,  0,  -  mg) R 


(21a) 


where  the  gravitational  acceleration  g  is  constant  for  range  firings. 
Substituting  the  right-hand  side  of  (21a)  in  (20),  we  see  that  the 
fixed-plane  components  of  C  depend  on  the  Euler  angle  0: 

£  =  (mg  sin  0,  0,  -  mg  cos  0)Rp  (21b) 


Note  that  by  the  definitions  of  the  two  Euler  angles,  the  angular 
velocity  of  the  intermediate  R'  (X'Y'Z')  system  with  respect  to  the 
range  system  is 


V(R)  =  (°* 
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Substituting  this  vector  in  (19) ,  we  obtain  the  angular  velocity  of  the 
fixed-plane  system  with  respect  to  the  range  system 


UFP(R) 


(-  $  sin  6,  6,  cos  6) 


FP 


(22) 


The  above  discussion  gives  us  the  physical  interpretation  of  \p  and 
6,  but  in  f'rder  to  work  with  these  angles,  we  needed  explicit  equations 
relating  \p  and  6  to  the  given  yaw  components  5  and  £  .  By  some  elemen- 

2  3 

tary  but  cumbersome  vector  analysis,  it  can  be  shown  that  the  desired 
relations  are 


where  i(»M  and  0M  are  the  missile's  pitch  and  yaw  angles,  respectively 
(see  Figure  2b)  : 

*m  =  sin'1  (-  v*)  *  sin_I  (y 

0  =  tan*1  (  — 3-  j  =  -  sin'1  |  - 3__  | 

M  \UJ  V005^/ 

where  aR  and  ER  are  the  azimuth  and  elevation  angles,  respectively,  of 
the  velocity  vector  in  the  range  system  (see  Figure  2c) : 


•  •  • 

where  x,y,z  are  the  velocity  components  in  the  range  system. 


(25) 

(26) 
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For  the  present  analysis  of  range  firings,  we  made  the  simplifying 
assumption  that  the  angles  aR  and  ER  could  be  ignored.  Letting 

aR  =  ER  =  0,  equations  (23)  and  (24)  reduce  to 

^  =  if>M  (23a) 

0  =  0M  (24a) 


By  equations  (17),  (25)  and  (26),  the  magnitude  6  of  the  yaw  can  be 
expressed  in  terms  of  the  pitch  and  yaw  angles: 


=  (sin2  i/»M  +  cos2 


♦m  sin“ 


ej 


(29) 


IV.  THE  EQUATIONS  OF  MOTION 

In  this  section,  the  working  forms  of  the  equations  of  motion  are 
derived  (albeit  briefly)  from  the  basic  expressions  for  Newton's  Second 
Law.  By  writing  out  this  derivation,  we  can  point  out  where  and  what 
kinds  of  assumptions  and  simplifications  were  made  and  thus  facilitate 
future  changes. 

The  interdependence  of  the  various  equations  and  of  the  three 
distinct  reductions  (drag,  yaw  and  CG)  are  indicated  in  Figure  3. 

A.  The  Drag  Equation 

The  classic  drag  equation  has  the  form 


where 


[MLT'2] 


■  drag  force 
=  -  (m  Aj  CD  V)  V 

\  ■  Hr1  i^i 

2  =  gravitational  force 


(30) 
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5  =  position  vector  of  the  missile's  CG 


and  where  the  subscript  I  denotes  vector  differentiation  in  an  inertial 
system.  For  enclosed-range  studies,  the  density  p  and  hence  A^  are 

known  constants.  The  drag  coefficient  C^  is  in  general  a  function  of 

Mach  number  and  c:  62.  For  each  of  the  twelve  Transonic  Range  rounds 
studied  here,  however,  we  could  ignore  the  Mach  number  variation  over 
the  observed  trajectory  and  assume  a  linear*  dependence  of  on  62: 


C 


D 


(31) 


We  seek  the  X-component  of  (30)  in  the  range  system.  For  the  pur¬ 
pose  of  performing  the  drag  reduction,  we  can  assume  that  the  range 
system  is  an  inertial  system  and  hence  ignore  the  centripetal  and 
Coriolis  accelerations  that  arise  in  any  earth-fixed  system.  Then  the 


X-component  of  (30)  can  be  written  as 

,, 

• 

~\ 

X 

II 

1 

> 

C„Vx 

1 

D 

=  -  A 

Cnv  V2 

► 

1 

DX 

where 

CDX  =  (?)  CD  =  S  COS  ER  C0S  aR 
=  down-range  drag  coefficient 

For  the  present  Transonic  Range  studies,  (32)  has  the  disadvantage 
that  the  independent  variable  t  -  which  should  be  known  exactly  at  each 
point  -  is  obtained  only  at  about  two-thirds  of  the  spark  stations  (and 
obtained  with  sufficient  accuracy  at  a  considerably  smaller  fraction) . 
On  the  other  hand,  the  down-range  coordinate  x  of  the  missile's  CG  is 
usually  known  very  accurately  at  each  station.  Thus  a  reasonable 
course  of  action  is  to  convert  the  independent  variable  in  (32)  from  t 
to  x.  We  have 


'Provision  was  made  in  the  coding  to  hca  die  C D  as  a  quadratic  in  f>2, 
but  the  higher-order  term  was  never  needed. 
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The  presence  in  (33)  of  the  velocity  V  is  a  nuisance  for  as  yet  we  have 
no  equation  for  generating  V.  For  the  nearly  horizontal  flights  en¬ 
countered  in  ballistic  ranges,  the  distinction  between  V  and 

x  (=  V  cos  E  cos  a  )  and  between  C  and  Cnv  can  be  ignored.  Thus  we 
K  K  L)  L)a 

can  write  (33)  in  the  final  form: 


Fhe  drag  reduction  -  that  is,  the  application  of  the  Chapman-Kirk 
technique  to  (34)  -  is  normally  done  before  the  yaw  reduction  and  so  at 
this  stage  we  know  the  values  of  62  only  at  certain  discrete  points 
(namely,  at  each  spark  station).  From  these  scattered  values,  we  ob¬ 
tained  additional  input  values  of  6 2  at  selected  values  of  x  by  inter¬ 
polation.  (After  the  yaw  reduction  has  been  performed,  the  drag 
reduction  could  be  re-done,  using  the  fitted  values  of  62  rather  than 
the  raw  plus  interpolated-raw  values.  For  the  present  study,  this 
wasn't  necessary.) 

In  addition  to  the  yaw  data,  the  required  input  for  a  given  round 
consisted  of  the  measured  times  t,  the  measured  x  values  and  initial 
estimates  of  the  two  explicit  parameters  j  and  the  two  initial 

conditions  (t  ,  V  ),  where  we  defined  initial  conditions  as  the  con- 
o  o 

ditions  at  the  first  spark  station.  The  output  consisted  of 

a.  the  least  squares  values  of  the  explicit  parameters  and  initial 
conditions ; 

b.  appropriate  error  estimates  in  the  above; 

c.  a  reliable  correlation  of  time  with  distance,  so  that  in  the 
remaining  equations  of  motion,  time  could  be  used  as  the  independent 
variable.  Of  course,  it  was  not  absolutely  necessary  to  revert  to  a 
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time  base.  For  theoretical  work,  the  arclength  r  along  the  trajectory 
°f  the  missile's  CC  (or  a  nondimensional  length  r/d)  is  often  a  more 
convenient  independent  variable  than  either  time  or  down-range  distance. 
The  equations  of  motion  (in  particular,  the  yaw  equations)  assume 

e 

simpler  forms  when  r  is  the  independent  variable.  Since  r  =  V,  our 
• 

approximation  x  a  V  is  equivalent  to  ignoring  the  distinction  between 
x  and  r.  We  could,  then,  have  written  all  the  equations  with  respect 
to  this  convenient  length  variable  which  we  measure  as  x  but  are  free 
to  interpret  as  r.  One  reason  we  didn't  is  that  the  present  study, 
while  interesting  in  itself,  is  also  regarded  as  preliminary,  getting- 
our-feet-wet  training  for  more  ambitious  studies  that  will  require  a 
time-based  set  of  equations.  By  working  with  time-based  equations  in 
the  present  study,  we  could  save  considerable  coding  effort.  Compre¬ 
hensive  computer  programs,  capable  of  handling  the  larger  problems, 
were  able  to  handle  the  current  problem  as  a  special  case. 

One  additional  bit  of  information  can  be  gleaned  from  our  reduction: 
a  representative  Cp  value  for  each  round,  obtained  by  replacing  62  in 

(31)  by  its  mean  value: 


C 


D 


P 


(35) 


Although  a  value  for  P  was  available  for  each  round  from  the  BRL  yaw 

reduction,  the  quantities  <*j.  =  sin'^y  ^jlisted  in  Tables  II  -  VI  were 

obtained  by  a  slightly  different  averaging  process  than  used  by  BRL. 
Thus  these  listed  values  differ  slightly  (by  less  than  4%)  from  the 
BRL  values.  Each  representative  drag  coefficient  CQ  can  be  compared 

with  the  range  value,  C  D  obtained  at  BRL  by  fitting  time  as  a  cubic 

in  down-range  distance  x: 


t 


C 


DR 


-'v 

t*  +  a  (x  -  x*)  +  a  (x  -  x*)2  +  a  (x  -  x*)3 
1  3 

2  a 

_ 2_ 

A  a 
1  1 


range  value  of  Cp  .it  [<  t,  x,  V) 


(f 


r  (36) 


J 


where  x*  is  the  given  mid-range  value  of  x. 


As  we  will  see,  the  variabie  V  appears  in  the  yaw  and  CG  equations 
and  the  question  arises:  how  should  we  generate  it  there?  We  could. 
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of  course,  store  the  large  number  of  discrete  velocity  values  available 
when  we  integrate  the  drag  equation  (34)  numerically.  This  would  be 
both  tedious  and  unnecessary.  Instead,  we  made  a  reasonable  assumption: 
for  purposes  of  performing  the  yaw  reduction*  on  the  given  Transonic 
Range  rounds,  V  can  be  represented  adequately  over  the  observed  tra¬ 
jectory  by  a  known  quadratic  in  time: 


V  +  V  (t  -  t  )  ♦  V  (t  -  t  )■ 
o  i  v  o'  2  0 


(37) 


Values  for  V  and  V  in  (37)  can  be  obtained  in  various  ways.  For 
1  2 

example,  if  we  replace  Cp  with  Cp  in  the  drag  equation,  the  solution 
can  be  written  at  once: 


V  =  VQ  exp  [-  Ai  CD  (x  -  x0)] 


1  +  CA  (t  ‘  V 


where 

““A  ’  "j  WD  *o 
with  the  truncated  series  expansion: 


C*  =  A,  Cn  V  [T'l] 


(38) 


V  =  V0  [1  -  CA  (t  -  t0)  *  Cl  (t  -  to)2]  (37a) 


B.  The  Roll  Equation 

Consider  a  missile-fixed  coordinate  system  XXX,  where  *he  X  - 

14  5  1 

axis  lies  along  the  missile's  longitudinal  axis  (as  in  the  XXX  fixed- 

12  3 

plane  system)  and  where  the  X  and  X  axes  are  rigidly  attached  to  the 

missile  in  a  right-handed  system.  Then  the  angular  velocity  of  this 
missile-fi/ed  system  relative  to  the  fixed-plane  system  is  given  by 


“MF(FP)  =  °’  °^FP 


(39) 


*The  assumption  is  not  necessary  for  the  CG  equations  where,  as  we  shall 
see,  V  can  be  generated  handily. 
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where  4>  is  the  roll  angle  (the  angle  between  the  X^X^  axes  and  the 

axes) .  By  (22)  and  (39) ,  we  can  write  an  expression  for  the  angular 
velocity  of  the  missile-fixed  system  with  respect  to  the  range  system: 


“MFCR)  =  “mF(FP)  +  WFP(R) 


=  (P,  0,  COS  0)Fp 


r  (40) 


where 


P 


-  tp  sin  9 


(41) 


The  angular  momentum  of  a  missile  with  rotational  symmetry  is  then 
t  e  (L  ,  L  ,  L  )pp  =  (Ix  p,  Iy  6,  Iy  *  cos  0)Fp  [ML2?'1]  (42) 


The  sum  of  the  moments  acting  on  the  missile  is  equal  to  the  time  deriv¬ 
ative  of  the  angular  momentum: 


M  =  (M  ,  M  ,  M  ) 


dt 


1*  2’  3JFP  “  dt 


•  •  • 


[ml2t‘2] 


-  [  (L^ ,  L^,  L^)  +  upp^Rj  X  ^]pp  W43) 


where  again  the  range  system  is  considered  an  inertial  system.  Hence 
by  (22)  and  (42), 


M  =  I  p 
1  x  ** 


(44) 


M  =  I0+(Ip+Iip  sin  0)  ip  cos  0 

2  y  *  y 


(45) 


M  =  I  i)i  cos  6-(I  p  ♦  2  I  iji  sin  0)  0 
3  y  x  v  y  y 


(46) 


Equation  (44)  is  the  roll  equation;  (45)  and  (46)  constitute  the  yaw 
equations . 

The  roll  equation  does  not  depend  on  the  other  two  moment  equations 
and  hence  can  be  solved  separately.  We  define  the  axial  moment  as 
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Mi  "  (p  t)  (t-)  (d)  )  Cip  [ml2j'Z]  (47) 


where  C 

X 

is  the  roll 

damping 

moment  coefficient*. 

Then 

the  roll  equa 

P 

tion  (44) 

becomes 

• 

P 

-  A  p  V 

3 

[T-2] 

(48) 

where 

A 

3 

■  W )  % 

[L-l] 

For  the  present  study,  A  can  be  considered  constant.  If  we  were  given 

3 

sufficient  spin  data  (obtained,  say,  by  measuring  the  position  on  each 
spark  photograph  of  two  distinguishable  pins  placed  in  the  base  of  the 
missile),  we  could  obtain  the  "best"  values  of  pQ  and  by  the  Chapman- 

Kirk  technique  or  by  a  fit  of  the  data  to  the  solution  of  (48) : 

P  =  P0  exP  [a3  (x  -  XD)]  (49) 

In  the  present  study,  however,  such  spin  data  was  unavailable.  Yet  the 
roll  rate  p  was  needed  in  the  yaw  '•quations.  We  resolved  this  problem 
by  assuming  that  p  is  a  known  linear  function  of  time: 

p  =  P0  [1  ♦  a3  Vo  (t  -  tQ)]  (50) 

where  Po  =  2  n  VQ/(nd),  rad/sec 

n  =  rifling,  cal/rev  (see  Table  I) 

and  where  A  was  evaluated  by  assuming  =  -  .013  for  all  twelve 
3  P 

rounds . 


*In  some  texts  (e.g.t  in  References  13  and  14), 
is  preferred  to  pd/V  and  in  those  texts  the 

P 


accordingly : 


(%  ■  t-) 


£d 

V 


the  combination  pd/(2V) 
must  be  interpreted 
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F 


C.  The  Yaw  Equations 

The  yaw  equations  describe  the  wobbling  motion  of  the  missile's 
longitudinal  axis  A  ^bout  the  tangent  to  the  trajectory,  that  i^,  about 
the  velocity  vector  V.  Any  two  variables  sufficient  to  orient  V  with 
respect  to  A  can  serve  as  the  dependent  variables.  In  the  usual  BRL 
Transonic  Range  reduction,  the  dependent  variables  are  the  yaw  components 
£  and  £  ;  here  we  use  the  Euler  angles  ip  and  9  and  the  governing  equa¬ 


tions  (45)  and  (46).  We  define  the  cross-moments  acting  on  the  missile 
as  follows*: 


^Provision  was  made  in  the  computer  program  to  cope  with  an  asymmetrical 
missile  by  including  some  trim  terms  not  shown  in  (51)  above.  These 
terns  were  not  used  in  the  present  study. 


**Note  that  C  in  (51)  is  multiplied  by  the  nondimensional  spin  pd/V, 
Mpa 

so  that  the  remarks  on  C  in  a  previous  footnote  apply  to  as  well. 


That  is. 


pa 


Likewise,  as  defined  here  is  half  the  of  References  13  aid  14. 

q  q 
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Substituting  (SI)  in  (45)  and  (46),  we  obtain  the  final  form  of  the  yaw 
equations : 


6 


<v  V  Si  + 

3  O 

- 


(U  p  d)  Cn  Md  V  8)  ^ 


pa 


[(Ix/Iy)  p  +  $  sin  e]  cos  6 


(53) 


*  = 


(V  u  )  C^j  +  (u  p  d)  CM  ♦  (d  V  Ip  cos  8) 


+  [(I  /I  )  p  +  2  iji  sin  6]  9  ^  /cos  6 
*  X 


pa 

0^-  /CO 


q 

(54) 


where  A  is  a  known  constant: 
2 


TT  p 


d3 


8  I 


[L-2] 


Equations  (53)  and  (54)  involve  four  initial  conditions,  the  eight  aero¬ 
dynamic  parameters  indicated  in  (52)  and  (possibly)  one  physical  param¬ 
eter,  I  /I  .  Thus  a  maximum  of  thirteen  constants  could  be  determined 
x  y 

for  each  round  by  an  application  of  the  Chapman-Kirk  technique  to  (53) 
and  (54).  However,  for  any  computer  run,  any  one  or  more  of  the  eight 
aerodynamic  constants  in  (52)  could  be  treated  as  a  known  constant 
rather  than  as  a  parameter  to  be  determined.  For  example,  was 

R2 

fixed  at  zero  (that  is,  ignored)  in  all  the  present  computer  runs. 

Note  that  in  addition  to  0,  ^  and  their  derivatives,  (53)  and  (54) 
involve  five  dependent  variables: 

V,  p,  62,  u  and  u 

2  3 

We  have  already  given  equations  adequate  for  simulating  V  and  p  as 
functions  of  time,  namely,  (37)  and  (50),  respectively.  The  squared 
yaw  can  be  obtained  by  (17) : 


(u2  +  u2)  /V2 
2  3 


(17a) 


provided  we  know  u  and  u  .  Thus  the  ability  to  solve  (53)  and  (54) 
2  3 

hinges  now  on  our  ability  to  generate  u  and  u  .  We  proceed  now  to 

2  3 

derive  the  required  equations.  Let 
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?  2  <".•  V  F3>FP 


=  the  aerodynamic  force  vector 


and  assume  that  the  only  forces  acting  on  the  missile  are  the  aerodynamic 
and  gravitational  forces.  Then  the  force  equation  is 


t  *  2  = 


d$ 

at 


•  • 


m  [Oy  u2.  u3)  ♦  -Fp(R)  x  VlFp 


or,  by  (21b)  and  (22) , 


(55) 


F  +  mg  sin  9  =  m  (u  -  u  cos  6  +  u  9) 

1  12  3 


(56) 


m  (u  +  u  i L  sin  8  u  \ji  cos  9) 
2  3  l 


(57) 


F  -  mg  cos  9  =  m(u  -  u  9-u  ^  sin  9) 

3  3  12 


(58) 


We  assume  that  the  fixed-plane  components  of  the  aerodynamic  force  are 
given  by 


F  =  -  m  A  C.  V  u 

1  1  D  1 


F  +  i  F 
2  3 


=  •  m  A,  (  S  V  +  1  CN  Pd  )  (u 
1  v  a  pa  / 


+  i  u  ) 

2  3 


V(59) 


where  for  the  present  purposes ,  and  Cw  are  considered  known 


N 


a  pa 

constants*.  Then  equations  (57)  ard  (58)  can  be  written  as: 


*As  tilth  , 
P 


£d 

V 
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where 


l.  [<V  V  \  -  (M  V  %.] 


-  (u  cos  6  +  u  sin  0) 
1  3 


u  =  -  A 


.  [<V  V  CN„  *  (P  d  V  V] 


•  • 


+  u  0  +  u  i p  sin  0  -  g  cos  0 
1  2 


(V2  -  u2  -  u2)5*  =  (1  -  fi2)  V 

2  3 


(60) 


(61) 


Thus  we  generate  u  and  u  by  assigning  values  to  C..  ,  C.,  ,  u  and  u 

2  3  Na  Npa  20  30 

and  then  solving  (60)  and  (61)  simultaneously  with  the  yaw  equations  (53) 
and  (54) . 

An  alternative  method  for  generating  u  and  u  follows  at  once  from 
assumptions  (23a)  and  (24a)  that  p  =  and  6  =  0M>  By  (16),  (25)  and 
(26)  we  have 


u  =  -  V  sin  i/< 

2 


u  =  V  cos  t|i  sin  0 
3 


These  expressions  eliminate  the  need  for  (60)  and  (61) .  The  only  draw¬ 
back  is  that  the  resulting  computer  program  couldn't  be  used  in  future 
cases  where  assumptions  (23a)  and  (24a)  are  invalid. 

While  we  didn't  use  this  simplification,  we  did  ease  the  labor  of 
computation  by  modifying  the  auxiliary  set  of  yaw  equations  produced 
by  parametric  differentiation"!  That  is,  we  replaced  the  functions  f ^ 

in  equation  (14)  by  approximating  expressions .  The  validity  of  this 
procedure  depends,  of  course,  on  the  adequacy  of  the  approximations. 

One  final  remark:  a  by-product  of  the  yaw  reduction  is  6 2  as  a 
function  of  time.  As  mentioned  in  section  IV  A,  a  second  drag  reduc¬ 
tion  could  be  performed  with  this  computed  62  as  input.  If  the  out¬ 
put  of  the  second  drag  reduction  differs  significantly  from  the 
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original  drag  output,  a  second  yaw  reduction  based  on  the  new  drag  data 
should  then  be  run.  For  the  present  study,  none  of  this  was  necessary. 

D.  The  CG  Equations 

The  motion  of  the  missile's  CG  along  the  trajectory  is  defined  by 
the  force  equation: 


?  *  t 


(62) 


We  are  interested  in  the  range  components  of  (62) .  Whereas  in  (30)  and 
(55)  wo  could  treat  the  range  system  as  an  inertial  system,  here  such 
simplification  is  not  as  defensible.  For  the  range  (earth-fixed) 
system,  the  acceleration  is  given  by 


_dM 

dt2 


I 


[(x,y,z)  +  ?  +  tj.]R 


(63) 


where  is  the  centripetal  acceleration,  which  we  promptly  ignore,  and 
where  2  is  the  Coriolis  acceleration,  which  we  retain: 

£  =  2  x  ^ 

where 

Ug  =  the  earth's  angular  velocity  vector 

=  wE  (cos  cos  9A,  cos  0^  sin  ©A,  sin  ©L)R 
=  2  it  radians/sidereal  day 

as  0.00007292  rad/sec 

and  where  @A  and  0^  are  known  constants  defined  in  the  List  of  Symbols. 
If  we  make  the  approximation 

$  =  (x,y,z)R  as  (V,  0,  0)R 
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then  the  Coriolis  acceleration  is  given  by 

£  =  2  «E  V  (0,  sin  0L>  -  cos  sin  0A)R  (64b) 

Defining 

?  =  (Fx,  Fy ,  Fz)r  =  CFj ,  Pa.  F^pp 


and  using  (21a),  (63)  and  (64b),  we  can  write  (62)  in  the  form 
x  =  Fx/m 

y  =  Fy/m  -  (2  wE  sin  0^)  V 
z  =  F^/m  +  (2  cos  0^  sin  0A)  V  -  g 


K65) 


To  obtain  expressions  i^,  Fy  and  F^,  we  assume  that  the  fixed-plane 
components  of  ?  have  the  form  already  indicated  in  (59) : 

F  =  -mA  CnVu 
1  1  D  ! 


F  +  i  F 
2  3 


m  A  CnA  V2 
l  DA 


-  m  A 


x  (  S  v  +  1  S  pd) 

1  \  a  pa  / 


(u  +  i  u  ) 

2  3 


r(66) 


where 


=  axial  drag  coefficient 


For  the  yaw  reduction,  C„  and  C.,  were  considered  known  constants  and 

N  N 

a  pa 

F  wasn't  needed;  here  we  assume 
1 
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CDA  =  CDAo  +  CDA2  6< 


C,  =  C,  +  C,  6' 


N  N 


ao  a  2 


=»  ■  S<  *S<  42 

pa  pao  pa2 


where  the  six  coefficients  on  the  right-hand  side  are  unknown  constants 
to  be  determined  by  the  Chapman-Kirk  technique*.  The  range  components 
of  the  aerodynamic  force  are  obtained  by  multiplying  the  fixed-plane 
components  by  the  transpose  of  the  transformation  matrix  given  in 
equation  (20) : 

F„  =  F  cos  6  cos  ip  -  F  sin  ip  *  F  sin  9  cos  <b  1 
X  1  2  3  i 


F 


Y 


F  cos  9  sin  ^  +  F  cos  ip  +  F  sin  6  sin  ip 
1  2  3 


>*  (67) 


F_  =  -  F  sin  9  +  F  cos 

z  1  3 

Substituting  (66)  and  (67)  in  (65) 
equations : 


8  J 

we  obtain  the  final  form  of  the  CG 

9  cos  ij; 


-  [or  V  SJ 

*  [o'  -3>  »2) 


sin  <|r 


sin  0  cos  ip  > 


(68) 


•Since  CM  -  (1  -  «*)  CD  •  .  (c^  -  C^)  52  ■  i\  chere 

and  Cn  are  well-determined  from  the  drag  reduction,  we  could,  as  an 

alternative,  consider  a  known  function  of  6  ,  thus  reducing  to  four 

the  number  of  aerodynamic  coefficients  to  be  determined  by  the  CG 
equations . 
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fws*? 


-  A  •< 


cos  6  sin  \p 


( 2  c“ ) 

'  [(V  .,)  ',B  - 

*  [lV  V  \  *  (>»*  V  V]  51 


(Pd  U  )  CN  I  C0S  * 

3  pa . 


1 


sin  0  sin  ij> 


-  (2  ujg  sin  0L)  V 


A  -< 


sin  8 


(v2  Sa)  sil 

(V  u  )  ♦  (pd  u  ) 

J  a  z  pa 

♦  (2  w  cos  0L  sin  0A)  V  -  g 


d 


cos  e 


) 


(69) 


(70) 


The  right-hand  sides  of  equations  (68  -  70)  involve  six  dependent 
variables : 


V,  p,  u^,  u^,  8  and  i|i. 

The  first  of  these  can  be  generated  internally: 

V  =  (x2  +  y2  +  z2)*5 

The  other  five  must  be  brought  in  from  outside.  For  the  spin  p,  we 

used  the  linear  expression  given  in  (50).  The  variables  u  ,  u  ,  8  and 

2  3 

i|>  were  obtained  from  the  yaw  reduction,  as  described  in  the  previous 
subsection. 

The  aerodynamic  coefficients  and  initial  conditions  of  equations 
(68  -  70)  were  adjusted  by  the  Chapman-Kirk  technique  until  the  solution 
of  the  equations  was  a  least  squares  fit  to  the  measured  (x,y,z)R  values 

of  the  CG. 
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V.  RESULTS 

The  results  obtained  in  the  present  study  can  be  interpreted  better 
in  the  light  of  a  little  background  on  the  four  shell  types. 

A.  Background 

The  90mm  M71  HE  shell  was  designed  to  be  fired  from  a  1/32  twist 
gun  at  Mach  2.4;  at  lower  muzzle  velocities,  the  shell  was  known  to 
perform  unsatisfactorily.  The  two  rounds  studied  here  (round  2-4203 
at  7°  average  yaw  and  round  2-4204  at  5°)  were  fired  at  about  Mach  0.93. 
At  that  muzzle  velocity  and  a  1/32  twist,  the  gyroscopic  stability  factor 
s^  was  slightly  less  than  unity  at  launch,  so  that  the  shell  was  gyro- 

scopically  unstable  when  it  emerged  from  the  tube.  However,  the  situ¬ 
ation  (like  the  shell)  took  some  turns  for  the  better.  As  the  shell 
travelled  downrange,  (a)  the  velocity  decay  exceeded  the  spin  decay  and 
(b)  the  static  moment  coefficient  decreased  with  decreasing  Mach  number. 
As  a  result,  s^  increased,  eventually  crossing  the  value  1.0,  so  that 

the  shell  became  and  remained  gyros copi cal ly  stable. 

The  M71  rounds  were  fired  and  initially  reduced  in  1956,  but  the 
values  obtained  for  some  of  the  coefficients  were  known  to  be  nonsense. 
The  photographic  plates  were  reread  several  times  in  the  next  four  years 
in  a  largely  unsuccessful  attempt  to  improve  the  linear  theory  fit.  One 
positive  result  of  all  these  readings  is  that  the  probable  errors  of  the 
fit  in  the  present  study  are  smaller  for  the  two  M71  rounds  than  for  any 
of  the  other  rounds. 

Before  applying  the  Chapman-Kirk  technique  to  the  two  M71  rounds, 
we  simulated  their  motion  on  an  analog  computer,  using  partially  linear¬ 
ized  equations  of  motion.  This  side  investigation  served  two  purposes: 
it  satisfied  our  curiosity  as  to  the  behavior  of  a  shell  flying  along 
the  borderline  of  gyroscopic  stability  -  instability  ard  it  provided 
accurate  estimates  of  the  initial  conditions  required  by  the  Chapman- 
Kirk  technique.  For  the  other  three  shell  types,  howiver,  we  were 
satisfied  to  use  a  simple  digital  subprogram  to  obtain  the  needed  first 
estimates . 

The  M329A1  shell  without  extension  was  a  4.2-inch  (107mm)  spin- 
stabilized  mortar  shell  to  which  a  subcaliber  cylindrica  after-body  - 
a  boom  -  had  been  attached  at  the  base  (see  Figure  1).  This  boom, 
about  0.35  caliber  in  diameter  and  0.7  caliber  long,  was  used  to  carry 
the  ignition  and  propelling  charge  and  to  provide  needed  volume.  The 
boom  served  no  aerodynamic  purpose  and  it  was  hoped  at  the  time  that 
the  boom  would  have  negligible  aerodynamic  effect  on  the  shell's  per¬ 
formance.  Careful  experiments11  proved  otherwise. 
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The  M329AI  shell  with  extension  carried,  as  its  name  might  imply, 
an  extension  to  the  original  boom,  bringing  the  total  boom  length  to 
1.35  calibers  (see  Figure  1).  Tests11  comparing  the  M329A1  with  and 
without  the  extension  clearly  indicated  that  the  added  length  affected 
the  aerodynamics,  decreasing  the  stability  and  increasing  the  drag  and 
the  shell's  tendency  to  fly  erratically  now  and  then. 

The  M329A1E1  shell  differed  from  the  two  previous  types  in  that  it 
had  a  longer  ogive,  a  shorter  body,  a  boattail  rather  than  a  square 
base  and  a  shorter  over-all  length  (see  Figure  1).  The  boom  was  about 
0.31  caliber  in  diameter  and  0.75  caliber  long. 

B.  Drag  Results 

A  comparison  of  the  BRL  drag  results  with  those  of  the  GE  Chapman- 
Kirk  approach  is  given  in  Table  II.  The  last  column  lists  the  BRL 
values  of  and  for  each  of  the  foui  projectile  types,  obtained 

0  2 

by  a  least  squares  fit  of  the  data  for  each  round  of  the  given  type  to 
the  equation 


CDR  =  CD  +  CD  C^BRL  ^71') 

0  2 

Although  additional  rounds  of  each  oi  the  four  shell  types  were  avail¬ 
able  for  fitting  BRL  data  to  equation  ."1),  it  was  considered  a  fairer 
comparison  to  work  solely  with  the  twelve  given  rounds.  As  a  result, 
the  fits  for  the  M71  and  M329A1E1  -  where  there  are  only  two  rounds 
each  -  are  exact  and  unreliable.  To  make  matters  worse,  the  yaw  levels 
of  the  two  M71  rounds  are  roughly  the  same,  so  that  we  can  expect  the 
slope  to  be  poorly  determined.  In  fact,  only  for  type  M329A1  with 
extension,  where  there  are  five  rounds  (albeit  only  three  yaw  levels), 
can  we  legard  the  results  of  the  fit  with  any  semblance  of  confidence. 

In  spite  of  all  this,  the  BRL  fitted  values  of  and  and  the  GE 

o  2 

Chapman-Kirk  values  are  not  wildly  dissimilar  -  they  are,  so  to  speak, 
in  the  same  ball  park  -  and  penaps  that  is  all  we  can  expect  from 
such  a  meager  number  of  rounds . 

One  encouraging  feature  of  the  GE  lata  is  that  the  values  of 

o 

and  Cp  obtained  are  approximately  the  same  for  each  round  of  a  given 

2 

type.  This  seems  a  necessary  if  not  sufficient  condition  for  trusting 
the  results.  For  each  round,  we  can  also  compare  Cp  (from  equation  35) 

with  the  BRL  range  value  CpR  shewn  in  the  next  to  last  column  of  Table 

II.  Here  we  find  that  the  agreement  is  quite  good  when  the  yaw  is  small 
ana  poorer  when  the  yaw  is  large. 
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c. 


Yaw  Results 


M 


The  yaw  results  are  shown  in  Tables  III  -  VI.  The  BRL  values  of 
CM  and  shown  for  each  round  are  -  for  the  small  yaw  rounds  - 

a  q  pa 

the  values  obtained  by  the  standard  BRL  epicycle  yaw  reduction1’2.  For 
the  large  yaw  rounds  (8615,  8618,  8621  and  8721) ,  the  BRL  reduction 
values  were  corrected  to  compensate  for  the  fact  that  certain  geometrical 
terms  in  the  yaw  reduction  -  terms  which  are  only  significant  for  large 
yew  angles  -  are  ignored  in  the  BRL  reduction. 


Again  it  should  be  noted  that  all  the  coefficients 


and  CM 


q  pa 

shown  in  these  tables,  both  BRL  and  GE  values,  are  defined  (see  equation 
51)  as  half  the  coefficients  designated  by  the  same  symbols  in,  for 
example,  References  9,  13  and  14.  This  disparity  is  unfortunate  but 
necessary  if  the  present  paper  is  to  be  consistent  with  other  BRL  reports. 


For  most  of  the  rounds,  more  than  one  GE  case  per  round  is  shown  in 
the  tables.  These  cases  differ  in  the  selection  of  which  parameters  were 
fixed  (indicated  by  an  asterisk)  and  which  were  allowed  to  seek  out  their 
optimum  values.  As  might  be  expected,  the  more  parameters  fitted,  the 
smaller  (with  a  few  exceptions)  was  the  probable  error  of  the  fit.  For 
three  rounds  (8713,  8716  and  8981),  GE  assumed  in  one  case  per  round 
that  all  the  moment  coefficients  were  constant.  This  affords  a  direct 
comparison  with  the  BRL  values. 


The  last  row  in  Tables  III  -  VI  lists  the  four  coefficients  obtain¬ 
ed  by  fitting  the  BRL  data  for  each  round  of  the  given  type  to  the  equa¬ 
tions 


where  the  constant  62  is  an  effective  squared  yaw,  obtained  as  a  by¬ 
product  of  the  BRL  yaw  reduction.  The  previous  remarks  on  the  fit  to 
Equation  (71)  apply  with  equal  force  to  Equation  (72)  :  the  results  of 
the  fit  are  not  too  trustworthy  but  are  least  suspect  for  the  M329A1 
with  extension  (Table  IV) . 


Round  8618  and  to  a  lesser  extent  round  8621  (Table  IV)  revealed 
a  very  strong  Magnus  moment  nonlinearity.  In  fact,  by  a  separate 
computer  curve  fit  program,  we  determined  that  the  expansion 
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Si  =  Si  +Si  62  +  Si  64 

pa  pao  pa2  pa4 

♦  66  +  6 0  (73) 

Pa6  Pa8 

was  needed  in  these  two  rounds  to  obtain  a  reasonable  fit  to  the  given 
data.  Accordingly,  our  main  program  was  modified  to  allow  for  the  last 
terms  in  (73)  -  terms  that  are  missing  in  equation  (52).  Unfortu¬ 

nately,  we  were  unable  to  determine  values  for  the  two  new  coefficients; 
the  process  diverged  on  every  attempt.  That  is,  it  was  not  possible 
with  the  present  technique  and/or  the  available  data  to  determine  more 
than  three  terms  of  the  Magnus  moment  coefficient  expansion. 

This  situation  should  be  investigated  further.  A  first  step  would 
be  to  analyze  by  the  present  technique  the  output  of  a  six-degrees-of- 
freedom  program  that  simulates  the  motion  of  a  projectile  with  a  highly 
nonlinear  Magnus  moment.  In  this  way,  the  quantity  and  the  quality  of 
the  experimental  data  needed  to  obtain  a  satisfactory  fit  for  this  type 
of  nonlinearity  could  be  determined. 

D.  CG  Results 


The  CG  equations  (68-70)  were  applied  only  to  the  two  M329A1E1 
rounds  (see  Table  VII).  The  values  obtained  for  the  normal  force 
coefficient  are  in  good  agreement  with  each  other  and  with  the  BRL 
values.  The  Magnus  force  coefficient  (which  is  less  by  that  ubiquitous 
factor  of  two  than  the  coefficient  labelled  Cy  in  Reference  14)  is 

pa 

not  very  well  determined.  The  BRL  values  of  the  axial  drag  coefficient 
shown  in  Table  VII  were  obtained  by  the  approximation 

Mm  ‘  C™ 

using  the  C  values  in  Table  III. 

UK 


VI.  CONCLUSIONS  AND  RECOMMENDATIONS 

We  have  shown  that  for  twelve  problem  rounds  the  Chapman-Kirk 
technique  could  be  applied  satisfactorily  to  free-flight,  enclosed- 
range  data.  However,  before  we  can  safely  apply  the  technique  on  a 
steady,  production  basis  to  high-yaw  rounds,  additional  theoretical 
study  is  needed.  A  systematic  investigation  should  be  made,  whereby  a 
variety  of  trajectories  are  computer-generated,  suitable  noise  is  intro¬ 
duced  and  the  resultant  data  fed  to  the  Chapman-Kirk  system. 
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Some  force  and  moment  coefficients  are  more  easily  determined  than 
others;  for  example,  those  depending  mainly  on  the  frequencies  of  the 
motion  are  more  easily  pinpointed  than  those  related  to  the  damping. 

For  a  given  force  or  moment  expansion  in  62,  the  coefficients  naturally 
become  less  reliable  as  the  order  of  the  term  increases  and  shortly  a 
limit  is  reached;  when  any  term  of  order  higher  than  this  limit  is  in¬ 
cluded  in  the  analysis,  the  process  fails  to  converge  (or,  worse  yet, 
converges  to  wrong  answers) .  This  limit  is  a  function  of  the  number  and 
accuracy  of  the  observations,  so  that  any  improvement  in  these  factors 
(up  to  some  point  of  diminishing  return)  would  be  helpful. 

It  might  be  possible  to  determine  high  order  terms  more  accurately 
for  a  given  set  of  data  if  the  lowest  order  terms  are  fixed  at  good 
values.  In  particular,  low  yaw  rounds  (say,  ci^  <  3°)  should  establish 

adequate  zero-yaw  coefficients  for  use  with  the  high  yaw  rounds.  In 
general,  the  fewer  the  number  of  coefficients  to  be  determined,  the 
less  likely  that  computational  instabilities  will  arise. 
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NOTE:  ALL  DIMENSIONS  ARE  IN  CALIBERS 


FIGURE  I  SIMPLIFIED  SKETCH  OF  THE  FOUR 
PROJECTILE  TYPES 
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FIG.  2o  f  IS  POSITIVE  FOR  A  CCW  ROTATION  OFTHE  XY- PLANE 
ABOUT  Z  (THUS  *  AS  SHOWN  IN  THE  FIGURE  IS  POSITIVE) 

0  IS  POSITIVE  FOR  A  CCW  ROTATION  OF  THE  Z'X'-PLANE 
ABOUT  Y1  (THUS  0  AS  SHOWN  IN  THE  FIGURE  IS  NEGATIVE) 
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Table  I.  Physical  Parameters  and  Environment 


Preceding  page  blank 


Table  II.  Dra^,  Reduction  Results 
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Table  I\ .  M329A1  (w/ext)  Yaw  Reduction  Results  (continued) 
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Table  V.  M329A1  (wo/ext)  Yaw  Reduction  Results 
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Table  VII.  M329A1E1  CG  Reduction  Results 
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